---
title: "Speciale"
author: "Jeppe Sabroe Thegen"
date: "4/11/2022"
output: pdf_document
editor_options: 
  chunk_output_type: console
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
options(scipen=999)

```

```{r Loads packages}
library(psych)
library(readxl)
library(tikzDevice)
library(cowplot)
library(tidyverse)
library(stargazer)
library(xtable)
library(missForest)
library(specr)
library(knitr)
library(nnet)
library(reshape2)
library(RColorBrewer)

```

```{r #Defining custom functions}
#Create function to normalize indexes between 0 and 1
range01 <- function(x, ...){(x - min(x, ...)) / (max(x, ...) - min(x, ...))}

#Create function to reverse variables
reverso <- function(x) {
y <- (x * (-1)) + max(x, na.rm = T)
return(y)
}

#Create a function to get the mode
getmode <- function(v) {
   uniqv <- unique(v)
   uniqv[which.max(tabulate(match(v, uniqv)))]
}
```

```{r Imports data}
setwd("~/Library/Mobile Documents/com~apple~CloudDocs/Speciale/Data")
drpanel <- read_excel("220009 - Klimaprojekt datatræk1.xlsx")
```

# Cleaning and munging

```{r Cleans independent variables}
#Makes gender a factor
drpanel$gender <- drpanel$`5_koen` %>% as.factor

#Makes region a factor
drpanel$region <- drpanel$`6_region` %>% as.factor

#Makes education an ordered factor
drpanel$edu <- factor(drpanel$`7_uddannelse`, ordered = TRUE, 
                                levels = c("Folkeskole 7 år eller derunder",
                                           "Folkeskole 8 - 10 år",
                                           "Erhvervseksamen",
                                           "Studentereksamen (inkl. HF, HH, HG, HTX)",
                                           "Kort videregående uddannelse (under 3 år)",
                                           "Mellem videregående uddannelse (3 til 4 år)",
                                           "Lang videregående uddannelse (5 år eller derover)"))

#Makes a numeric variable for education years
drpanel$eduyears <- drpanel$edu
drpanel$eduyears <- dplyr::recode(drpanel$eduyears, 
                                  "Folkeskole 7 år eller derunder" = 7,
                                  "Folkeskole 8 - 10 år" = 10,
                                  "Erhvervseksamen" = 12,
                                  "Studentereksamen (inkl. HF, HH, HG, HTX)" = 13,
                                  "Kort videregående uddannelse (under 3 år)" = 16,
                                  "Mellem videregående uddannelse (3 til 4 år)" = 17,
                                  "Lang videregående uddannelse (5 år eller derover)" = 18)

#Makes income an ordered factor
drpanel$income <- factor(drpanel$`Personlig Indkomst`, ordered = TRUE, 
                                levels = c("Under 100.000 kr.",
                                           "100.000 kr. – 299.999 kr.",
                                           "300.000 kr. – 499.999 kr.",
                                           "500.000 kr. – 699.999 kr.",
                                           "700.000 kr. eller derover"))



#Only 14 respondents have less than 7 years of primary school. I combine them with the rest, as is normal procedure in Danish surveys.

levels(drpanel$edu)[levels(drpanel$edu)=="Folkeskole 7 år eller derunder"] <-"Folkeskole 8 - 10 år"

#Makes party an un-ordered factor
drpanel$party <- 1:1518

for(i in 1:nrow(drpanel)){
  
    if(drpanel$Background3[i] %in% 12) {
    drpanel$party[i] <- drpanel$`Andet parti, notér:Background3`[i]
  } else {
    drpanel$party[i] <- drpanel$Background3[i]}
}


#Some want to vote for Moderaterne or the Vegan Party. Enough to make parties

moderatevoters <- grep('Moderat|løkke', drpanel$party, ignore.case=T)
drpanel$party[moderatevoters]<- "Moderaterne"

veganvoters <- grep('vegan', drpanel$party, ignore.case=T)
drpanel$party[veganvoters]<- "Veganerpartiet"

#Code as factor
drpanel$party[drpanel$party == 1] <- "Socialdemokratiet"
drpanel$party[drpanel$party == 2] <- "Radikale Venstre"
drpanel$party[drpanel$party == 3] <- "Det Konservative Folkeparti"
drpanel$party[drpanel$party == 4] <- "Nye Borgerlige"
drpanel$party[drpanel$party == 5] <- "Socialistisk Folkeparti"
drpanel$party[drpanel$party == 6] <- "Liberal Alliance"
drpanel$party[drpanel$party == 7] <- "Kristendemokraterne"
drpanel$party[drpanel$party == 8] <- "Dansk Folkeparti"
drpanel$party[drpanel$party == 9] <- "Venstre"
drpanel$party[drpanel$party == 10] <- "Enhedslisten"
drpanel$party[drpanel$party == 11] <- "Alternativet"
drpanel$party[drpanel$party %in% c(14,15)] <- "Stemmer ikke/Stemmer blankt"
drpanel$party[drpanel$party %in% c(13,16,98,99)] <- "Stemmer ikke/Stemmer blankt"

#Could do one where blanks are missing instead.

#Codes left/right at this point before loss of detail by dropping parties with only one vote.
drpanel$right <- 1:nrow(drpanel)
for(i in 1:nrow(drpanel)){
    if(drpanel$party[i] %in% c("Socialdemokratiet","Radikale Venstre","Socialistisk Folkeparti","Enhedslisten","Alternativet", "Veganerpartiet", "Grønlandske politikere på venstrefløjen", "Klimapartiet momentum", "Kommunisterne")) {
    drpanel$right[i] <- 0} 
  else if (drpanel$party[i] %in% "Stemmer ikke/Stemmer blankt") {
    drpanel$right[i] <- NA
  } else {
    drpanel$right[i] <- 1}
}

#Respondents are markedly more red (807/532) than the Danish population.

#Now we can parties with only one voter from the data.
drpanel$party[drpanel$party %in% c("De Grønne", "Klimapartiet momentum", "Grønlandske politikere på venstrefløjen", "Kommunisterne")] <- NA

#Codes ideology between 0 and 1. One is right wing
drpanel$ideology <- range01(drpanel$Background2, na.rm=T)


#Codes treatments in detail
drpanel$treatment <- factor(drpanel$Samplenøgle, ordered = FALSE, 
                                levels = c("Kontrol",
                                           "Geografisk",
                                           "Social",
                                           "Risiko",
                                           "Tid",
                                           "Kontrol+håb",
                                           "Geografisk+håb",
                                           "Social+håb",
                                           "Risiko+håb",
                                           "tid+håb",
                                           "Ingen Treatment"))

#Codes treatments more broadly
drpanel$treatment5 <- drpanel$treatment
levels(drpanel$treatment5)[levels(drpanel$treatment5)=="Risiko+håb"] <-"Risiko"
levels(drpanel$treatment5)[levels(drpanel$treatment5)=="Kontrol+håb"] <-"Kontrol"
levels(drpanel$treatment5)[levels(drpanel$treatment5)=="tid+håb"] <-"Tid"
levels(drpanel$treatment5)[levels(drpanel$treatment5)=="Geografisk+håb"] <-"Geografisk"
levels(drpanel$treatment5)[levels(drpanel$treatment5)=="Social+håb"] <-"Social"

#Codes treatments as hope or not
drpanel$treatment_hope <- drpanel$treatment
levels(drpanel$treatment_hope)[levels(drpanel$treatment_hope)=="Risiko+håb"] <-"Hope"
levels(drpanel$treatment_hope)[levels(drpanel$treatment_hope)=="Kontrol+håb"] <-"Hope"
levels(drpanel$treatment_hope)[levels(drpanel$treatment_hope)=="tid+håb"] <-"Hope"
levels(drpanel$treatment_hope)[levels(drpanel$treatment_hope)=="Geografisk+håb"] <-"Hope"
levels(drpanel$treatment_hope)[levels(drpanel$treatment_hope)=="Social+håb"] <-"Hope"

levels(drpanel$treatment_hope)[levels(drpanel$treatment_hope)=="Risiko"] <-"None"
levels(drpanel$treatment_hope)[levels(drpanel$treatment_hope)=="Kontrol"] <-"None"
levels(drpanel$treatment_hope)[levels(drpanel$treatment_hope)=="Tid"] <-"None"
levels(drpanel$treatment_hope)[levels(drpanel$treatment_hope)=="Geografisk"] <-"None"
levels(drpanel$treatment_hope)[levels(drpanel$treatment_hope)=="Social"] <-"None"

#Codes treatments as hope or not
drpanel$treated <- drpanel$treatment
levels(drpanel$treated)[levels(drpanel$treated)=="Risiko+håb"] <- 1
levels(drpanel$treated)[levels(drpanel$treated)=="Kontrol+håb"] <- 1
levels(drpanel$treated)[levels(drpanel$treated)=="tid+håb"] <- 1
levels(drpanel$treated)[levels(drpanel$treated)=="Geografisk+håb"] <-1
levels(drpanel$treated)[levels(drpanel$treated)=="Social+håb"] <-1

levels(drpanel$treated)[levels(drpanel$treated)=="Risiko"] <-1
levels(drpanel$treated)[levels(drpanel$treated)=="Kontrol"] <-1
levels(drpanel$treated)[levels(drpanel$treated)=="Tid"] <-1
levels(drpanel$treated)[levels(drpanel$treated)=="Geografisk"] <-1
levels(drpanel$treated)[levels(drpanel$treated)=="Social"] <-1

levels(drpanel$treated)[levels(drpanel$treated)=="Ingen Treatment"] <-0


drpanel$treated <- drpanel$treated %>% as.numeric

drpanel$treatment_hope_num <- drpanel$treatment_hope
drpanel$treatment_hope_num[drpanel$treatment_hope_num == "Ingen Treatment"] <- NA
drpanel$treatment_hope_num <- drpanel$treatment_hope_num %>% as.numeric()

#Just age in English for parallelism
drpanel$age <- drpanel$`4_alder`



```

```{r Cleans and prepares dependent variables as well as makes indexes}
drpanel[drpanel == 99] <- NA
drpanel[drpanel == 88] <- NA


#Concern 2 and 4 have to be reversed
drpanel$Concern2_reverse <- drpanel$Concern2 %>% reverso
drpanel$Concern4_reverse <- drpanel$concern4 %>% reverso

drpanel <- mutate(drpanel, concern_index = rowMeans(dplyr::select(drpanel,c("Concern1", "Concern2_reverse","concern3","Concern4_reverse")), na.rm = TRUE))

#Based on the factor analysis it seems these three strongly correlate
drpanel <- mutate(drpanel, concern_index3 = rowMeans(dplyr::select(drpanel,c("Concern1", "Concern2_reverse","concern3")), na.rm = TRUE))

#Normalizing
drpanel$concern_index3 <- range01(drpanel$concern_index3, na.rm = T)

#Normalizes variables concerning concern
drpanel$concern_index <- range01(drpanel$concern_index, na.rm = T)
drpanel$concern3 <- range01(drpanel$concern3, na.rm = T)
drpanel$Concern1 <- range01(drpanel$Concern1, na.rm = T)
drpanel$Concern2_reverse <- range01(drpanel$Concern2_reverse, na.rm = T)
drpanel$Concern4_reverse <- range01(drpanel$Concern4_reverse, na.rm = T)
drpanel$concern5 <- range01(drpanel$concern5, na.rm = T)

#Policy-support
drpanel <- mutate(drpanel, support_index = rowMeans(dplyr::select(drpanel, c("policysupport1","policysupport2","policysupport3")), na.rm = TRUE))

#Normalizing
drpanel$support_index <- range01(drpanel$support_index, na.rm = T)

#Intentions
drpanel$Intention3_reverse <- drpanel$Intention3 %>% reverso

drpanel <- mutate(drpanel, intention_index = rowMeans(dplyr::select(drpanel, c("Intention1","Intention2","Intention3_reverse", "intention4")), na.rm = TRUE))

#Normalizing
drpanel$intention_index <- range01(drpanel$intention_index, na.rm = T)

drpanel <- mutate(drpanel, intention_index3 = rowMeans(dplyr::select(drpanel, c("Intention1","Intention2", "intention4")), na.rm = TRUE))

#Hope
drpanel$hope2_reverse <- drpanel$hope2 %>% reverso

drpanel <- mutate(drpanel, hope_index = rowMeans(dplyr::select(drpanel, c("hope1","hope2_reverse","hope3")), na.rm = TRUE))

#Normalizing
drpanel$hope_index <- range01(drpanel$hope_index, na.rm = T)

drpanel <- mutate(drpanel, hope_index12 = rowMeans(dplyr::select(drpanel, c("hope1","hope2_reverse")), na.rm = TRUE))

#Normalizing
drpanel$hope_index12 <- range01(drpanel$hope_index12, na.rm = T)

drpanel <- mutate(drpanel, hope_index13 = rowMeans(dplyr::select(drpanel, c("hope1","hope3")), na.rm = TRUE))

#Normalizing
drpanel$hope_index13 <- range01(drpanel$hope_index13, na.rm = T)

#Efficacy
drpanel$efficacy3_reverse <- drpanel$efficacy3 %>% reverso

drpanel <- mutate(drpanel, efficacy_index = rowMeans(dplyr::select(drpanel, c("efficacy1","efficacy2","efficacy3_reverse", "efficacy4")), na.rm = TRUE))

#Based on the factor analysis, maybe leave item 3 out.
drpanel <- mutate(drpanel, efficacy_index3 = rowMeans(dplyr::select(drpanel, c("efficacy1","efficacy2", "efficacy4")), na.rm = TRUE))

#Normalizing
drpanel$efficacy_index <- range01(drpanel$efficacy_index, na.rm = T)
drpanel$efficacy_index3 <- range01(drpanel$efficacy_index3, na.rm = T)

#Scepsis
drpanel$sceptisism1_reverse <- drpanel$sceptisism1 %>% reverso

drpanel <- mutate(drpanel, scepticism_index = rowMeans(dplyr::select(drpanel, c("sceptisism1_reverse","sceptisism2","sceptisism3")), na.rm = TRUE))

#Normalizing
drpanel$scepticism_index <- range01(drpanel$scepticism_index, na.rm = T)

```


```{r Cleans the fear items and creates an index}
drpanel$bange <- range01(drpanel$`Bange:Fear`,na.rm=T)
drpanel$truet <- range01(drpanel$`Truet:Fear`,na.rm=T)
drpanel$nervøs <- range01(drpanel$`Nervøs:Fear`,na.rm=T)
drpanel$skræmt <- range01(drpanel$`Skræmt:Fear`,na.rm=T)
drpanel$bekymret <- range01(drpanel$`Bekymret:Fear`,na.rm=T)

drpanel <- mutate(drpanel, fear_index = rowMeans(dplyr::select(drpanel, c("bange","truet","nervøs", "skræmt","bekymret")), na.rm = T))

drpanel <- mutate(drpanel, fear_index4 = rowMeans(dplyr::select(drpanel, c("bange","truet","nervøs", "skræmt","bekymret")), na.rm = T))

drpanel$optimistisk <- range01(drpanel$`Optimistisk:Fear`,na.rm=T) %>% reverso
drpanel$håbefuld <- range01(drpanel$`Håbefuld:Fear`,na.rm=T)


```


```{r Shows distributions of indexes}
#For aesthetic

drpanel$dummy <- drpanel$`1_indid`
for (i in 1:nrow(drpanel)){
  drpanel$dummy[i] <- "var"
}

melted_concern <- melt(drpanel %>% dplyr::select(c("concern_index", "concern_index3")))

gg_concern_index <- ggplot(melted_concern, aes(x=value)) + geom_density(aes(color=variable), show.legend = F)+ geom_point(aes(y=value, color = variable), alpha=0) + geom_histogram(aes(y=..density.., fill=variable), alpha=0.3, bins = 20, position = 'identity') + geom_vline(aes(xintercept=mean(value[variable=="concern_index"], na.rm=T)), linetype = "dashed", color="#E41A1C") + 
  geom_vline(aes(xintercept=mean(value[variable=="concern_index3"], na.rm=T)), linetype = "dashed", color="#377EB8") +
  theme_minimal() + scale_fill_brewer(palette = "Set1",direction =) +  scale_colour_brewer(palette = "Set1", labels = c("Bekymring", "Bekymring*")) + theme(legend.position = "bottom", axis.title.x=element_blank()) + guides(color = guide_legend(override.aes = list(shape = c("|", "|"), alpha = 1, size = 4)),
           title = "none",
           fill = "none",
           linetype = 0) + labs(color="", y = "Tæthed")


#Need to show both
melted_efficacy <- melt(drpanel %>% dplyr::select(c("efficacy_index", "efficacy_index3")))

gg_efficacy_index <- ggplot(melted_efficacy, aes(x=value)) + geom_density(aes(colour=variable), show.legend=F) + geom_histogram(aes(y=..density.., fill=variable), alpha=0.3, bins = 19, position='identity') +
  geom_point(aes(y=value, color = variable), alpha=0) + geom_vline(aes(xintercept=mean(value[variable=="efficacy_index"], na.rm=T)), linetype = "dashed", color="#E41A1C") + geom_vline(aes(xintercept=mean(value[variable=="efficacy_index3"], na.rm=T)), linetype = "dashed", color="#377EB8") + theme_minimal() + scale_fill_brewer(palette = "Set1",direction =) +  scale_colour_brewer(palette = "Set1", labels = c("Efficacy", "Efficacy§")) + theme(legend.position = "bottom", axis.title.x=element_blank()) + guides(color = guide_legend(override.aes = list(shape = c("|", "|"), alpha =1, size = 4)),
           title = "none",
           fill = "none",
           linetype=0) + labs(color="", y = "Tæthed")


gg_fear_index <- ggplot(drpanel, aes(fear_index)) + geom_density(aes(color=dummy)) + geom_histogram(aes(y=..density..), alpha=0.5, bins = 26) + geom_vline(aes(xintercept=mean(fear_index, na.rm=T)), linetype = "dashed") +
  geom_point(aes(y=fear_index, color = dummy), alpha=0) + theme_minimal() + theme(legend.position = "bottom", axis.title.x=element_blank()) +
  guides(color = guide_legend(override.aes = list(linetype = 0,
                                                  alpha= 1, 
                                                  size = 4,
                                                  shape = "|"))) + 
  labs(color = "", y = "Tæthed") + scale_color_discrete(labels = "Frygt")

gg_intention_index <- ggplot(drpanel, aes(intention_index)) + geom_density(aes(color=dummy)) + geom_histogram(aes(y=..density..), alpha=0.5, bins = 26) + geom_vline(aes(xintercept=mean(intention_index, na.rm=T)), linetype = "dashed") + geom_point(aes(y=intention_index, color = dummy), alpha=0) + theme_minimal() + theme(legend.position = "bottom", axis.title.x=element_blank()) +
  guides(color = guide_legend(override.aes = list(linetype = 0,
                                                  alpha= 1, 
                                                  size = 4,
                                                  shape = "|"))) + 
  labs(color = "", y = "Tæthed") + scale_color_discrete(labels = "Intentioner")

melted_hope <- melt(drpanel %>% dplyr::select(c("hope_index", "hope_index12", "hope_index13")))

gg_hope_index <- ggplot(melted_hope, aes(value)) + geom_density(aes(color=variable), show.legend = F) + theme_minimal()  + geom_vline(aes(xintercept=mean(value[variable=="hope_index"], na.rm=T)), linetype = "dashed", color="#E41A1C") + geom_point(aes(y=value, color = variable), alpha=0) +
geom_vline(aes(xintercept=mean(value[variable=="hope_index12"], na.rm=T)), linetype = "dashed", color="#377EB8") + 
geom_vline(aes(xintercept=mean(value[variable=="hope_index13"], na.rm=T)), linetype = "dashed", color="#4DAF4A") + 
  theme_minimal() + scale_fill_brewer(palette = "Set1",direction =) +  scale_colour_brewer(palette = "Set1", labels = c("Håb", "Håb†", "Håb‡")) + theme(legend.position = "bottom", axis.title.x=element_blank()) + guides(color = guide_legend(override.aes = list(shape = c("|", "|", "|"), alpha = 1, size = 4)),
           title = "none",
           fill = "none",
           linetype = 0 ) + labs(color="", y = "Tæthed")

gg_support_index <- ggplot(drpanel, aes(support_index)) + geom_density(aes(color=dummy)) + geom_histogram(aes(y=..density..), alpha=0.5, bins = 19) + geom_vline(aes(xintercept=mean(support_index, na.rm=T)), linetype = "dashed") + geom_point(aes(y=scepticism_index, color = dummy), alpha=0) + theme_minimal() + theme(legend.position = "bottom", axis.title.x=element_blank()) +
  guides(color = guide_legend(override.aes = list(linetype = 0,
                                                  alpha= 1, 
                                                  size = 4,
                                                  shape = "|"))) + 
  labs(color = "", y = "Tæthed") + scale_color_discrete(labels = "Policy support")

gg_scepticism_index <- ggplot(drpanel, aes(scepticism_index)) + geom_density(aes(color=dummy), show.legend = T) + geom_histogram(aes(y=..density..), alpha=0.5, bins = 19, show.legend = F) + geom_vline(aes(xintercept=mean(scepticism_index, na.rm=T)), linetype = "dashed") +
  geom_point(aes(y=scepticism_index, color = dummy), alpha=0) + theme_minimal() + theme(legend.position = "bottom", axis.title.x=element_blank()) +
  guides(color = guide_legend(override.aes = list(linetype = 0,
                                                  alpha= 1, 
                                                  size = 4,
                                                  shape = "|"))) + 
  labs(color = "", y = "Tæthed") + scale_color_discrete(labels = "Skepticisme")


pg_dr_indekser <- plot_grid(gg_concern_index, gg_efficacy_index, gg_intention_index, gg_support_index, gg_hope_index, gg_fear_index, gg_scepticism_index, 
          ncol =2, align = "hv")

tikz(file = "DRindekser.tex",
     standAlone=F,
     width = 6,
     height = 8,
     sanitize =T,
     )

print(pg_dr_indekser)

dev.off()

```


```{r Cleans the manipulation checks}
drpanel$check_faraway <- range01(drpanel$`...klimaforandringerne sker langt væk:Check`, na.rm=T)

drpanel$check_likeme <- range01(drpanel$`...klimaforandringerne vil ramme folk som mig:Check`, na.rm=T) %>% reverso

drpanel$check_notlikely <- range01(drpanel$`...det er usandsynligt at klimaforandringer vil indtræffe:Check`, na.rm=T)

drpanel$check_longtime <- range01(drpanel$`...klimaforandringer vil indtræffe snart:Check`, na.rm=T) %>% reverso

drpanel$feels_hope <- range01(drpanel$`...der er håb for at gøre noget ved klimaforandringerne:Check`, na.rm=T)

```

# Factor analyses

```{r Runs a factor analysis of concern}
df_factor_concern <- dplyr::select(drpanel,c("Concern1","Concern2_reverse","concern3","Concern4_reverse", "concern5"))

df_factor_concern <- df_factor_concern[complete.cases(df_factor_concern), ]

factanal.concern <- factanal(df_factor_concern, factors=2, scores = c("regression"), rotation = "varimax", na.action = "na.omit")
print(factanal.concern)

#Scree
fafitfree <- fa(df_factor_concern,nfactors = ncol(df_factor_concern), rotate = "varimax")
n_factors <- length(fafitfree$e.values)
scree     <- data.frame(
  Factor_n =  as.factor(1:n_factors), 
  Eigenvalue = fafitfree$e.values)
screeplot_concern<- ggplot(scree, aes(x = Factor_n, y = Eigenvalue, group = 1)) + 
  geom_point() + geom_line() +
  xlab("Faktorer") +
  ylab("Eigenværdi") +
  labs(title = "Bekymring") +
  geom_hline(yintercept=1, linetype=3, size =1) +
  theme_minimal()+ theme(plot.title = element_text(hjust = 0.5))

#screeplot

```

```{r Runs a factor analysis of efficacy}
df_factor_efficacy <- dplyr::select(drpanel,c("efficacy1","efficacy2","efficacy3_reverse", "efficacy4"))

df_factor_efficacy <- df_factor_efficacy[complete.cases(df_factor_efficacy), ]

factanal.efficacy <- factanal(df_factor_efficacy, factors=1, scores = c("regression"), rotation = "varimax", na.action = "na.omit")
print(factanal.efficacy)

#Scree
fafitfree <- fa(df_factor_efficacy,nfactors = ncol(df_factor_efficacy), rotate = "varimax")
n_factors <- length(fafitfree$e.values)
scree     <- data.frame(
  Factor_n =  as.factor(1:n_factors), 
  Eigenvalue = fafitfree$e.values)
screeplot_efficacy<- ggplot(scree, aes(x = Factor_n, y = Eigenvalue, group = 1)) + 
  geom_point() + geom_line() +
  xlab("Faktorer") +
  ylab("Eigenværdi") +
  labs(title = "Efficacy") +
  geom_hline(yintercept=1, linetype=3, size =1) +
  theme_minimal()+ theme(plot.title = element_text(hjust = 0.5))

#screeplot

```

```{r Runs a factor analysis of policy support}
df_factor_polsup <- dplyr::select(drpanel,c("policysupport1","policysupport2","policysupport3"))

df_factor_polsup <- df_factor_polsup[complete.cases(df_factor_polsup), ]

factanal.polsup <- factanal(df_factor_polsup, factors=1, scores = c("regression"), rotation = "varimax", na.action = "na.omit")
print(factanal.polsup)
#Generelt høje loadings

#Scree
fafitfree <- fa(df_factor_polsup,nfactors = ncol(df_factor_polsup), rotate = "varimax")
n_factors <- length(fafitfree$e.values)
scree     <- data.frame(
  Factor_n =  as.factor(1:n_factors), 
  Eigenvalue = fafitfree$e.values)
screeplot_polsup<- ggplot(scree, aes(x = Factor_n, y = Eigenvalue, group = 1)) + 
  geom_point() + geom_line() +
  xlab("Faktorer") +
  ylab("Eigenværdi") +
  labs(title = "Støtter politisk handling") +
  geom_hline(yintercept=1, linetype=3, size =1) +
  theme_minimal()+ theme(plot.title = element_text(hjust = 0.5))

#screeplot

```

```{r Runs a factor analysis of intentions}
df_factor_intentions <- dplyr::select(drpanel,c("Intention1","Intention2","Intention3_reverse", "intention4"))

df_factor_intentions <- df_factor_intentions[complete.cases(df_factor_intentions), ]

factanal.intentions <- factanal(df_factor_intentions, factors=1, scores = c("regression"), rotation = "varimax", na.action = "na.omit")
print(factanal.intentions)
#Item3 somewhat low. I don't see any obvious reason to leave out.

#Scree
fafitfree <- fa(df_factor_intentions,nfactors = ncol(df_factor_intentions), rotate = "varimax")
n_factors <- length(fafitfree$e.values)
scree     <- data.frame(
  Factor_n =  as.factor(1:n_factors), 
  Eigenvalue = fafitfree$e.values)
screeplot_intentions<- ggplot(scree, aes(x = Factor_n, y = Eigenvalue, group = 1)) + 
  geom_point() + geom_line() +
  xlab("Faktorer") +
  ylab("Eigenværdi") +
  labs(title = "Handlingsintentioner") +
  geom_hline(yintercept=1, linetype=3, size =1) +
  theme_minimal()+ theme(plot.title = element_text(hjust = 0.5))

#screeplot

```

```{r Runs a factor analysis of fear}
df_factor_fear <- dplyr::select(drpanel,c("bange","truet","nervøs", "skræmt","bekymret"))

df_factor_fear <- df_factor_fear[complete.cases(df_factor_fear), ]

factanal.fear <- factanal(df_factor_fear, factors=2, scores = c("regression"), rotation = "varimax", na.action = "na.omit")
print(factanal.fear)

#Scree
fafitfree <- fa(df_factor_fear,nfactors = ncol(df_factor_fear), rotate = "varimax")
n_factors <- length(fafitfree$e.values)
scree     <- data.frame(
  Factor_n =  as.factor(1:n_factors), 
  Eigenvalue = fafitfree$e.values)
screeplot_fear<- ggplot(scree, aes(x = Factor_n, y = Eigenvalue, group = 1)) + 
  geom_point() + geom_line() +
  xlab("Faktorer") +
  ylab("Eigenværdi") +
  labs(title = "Frygt") +
  geom_hline(yintercept=1, linetype=3, size =1) +
  theme_minimal()+ theme(plot.title = element_text(hjust = 0.5))

#screeplot

```

```{r Runs a factor analysis of hope}
df_factor_hope <- dplyr::select(drpanel,c("hope1","hope2_reverse","hope3"))

df_factor_hope <- df_factor_hope[complete.cases(df_factor_hope), ]

factanal.hope <- factanal(df_factor_hope, factors=1, scores = c("regression"), rotation = "varimax", na.action = "na.omit")
factanal.hope


#Scree
fafitfree <- fa(df_factor_hope,nfactors = ncol(df_factor_hope), rotate = "varimax")
n_factors <- length(fafitfree$e.values)
scree     <- data.frame(
  Factor_n =  as.factor(1:n_factors), 
  Eigenvalue = fafitfree$e.values)
screeplot_hope<- ggplot(scree, aes(x = Factor_n, y = Eigenvalue, group = 1)) + 
  geom_point() + geom_line() +
  xlab("Faktorer") +
  ylab("Eigenværdi") +
  labs(title = "Håb") +
  geom_hline(yintercept=1, linetype=3, size =1) +
  theme_minimal()+ theme(plot.title = element_text(hjust = 0.5))

#screeplot

```

```{r Runs a factor analysis of scepsis}
df_factor_scepticism <- dplyr::select(drpanel,c("sceptisism1_reverse","sceptisism2","sceptisism3"))

df_factor_scepticism <- df_factor_scepticism[complete.cases(df_factor_scepticism), ]

factanal.sceptism <- factanal(df_factor_scepticism, factors=1, scores = c("regression"), rotation = "varimax", na.action = "na.omit")
factanal.sceptism


#Scree
fafitfree <- fa(df_factor_scepticism,nfactors = ncol(df_factor_scepticism), rotate = "varimax")
n_factors <- length(fafitfree$e.values)
scree     <- data.frame(
  Factor_n =  as.factor(1:n_factors), 
  Eigenvalue = fafitfree$e.values)
screeplot_scepticism<- ggplot(scree, aes(x = Factor_n, y = Eigenvalue, group = 1)) + 
  geom_point() + geom_line() +
  xlab("Faktorer") +
  ylab("Eigenværdi") +
  labs(title = "Skepsis") +
  geom_hline(yintercept=1, linetype=3, size =1) +
  theme_minimal() + theme(plot.title = element_text(hjust = 0.5))

#screeplot

```

```{r Saving those screes as Tikz for latex}

screegrid<- plot_grid(screeplot_concern, screeplot_efficacy, screeplot_polsup, screeplot_intentions, screeplot_fear, screeplot_hope, screeplot_scepticism, 
          ncol =2)

#ggdraw(screegrid) + draw_label("Draft", colour = "#80404080", size = 120, angle = 45, xlim = c(0, 1), ylim = c(0, 1))

# tikz(file = "my_output_file3.tex", 
#      standAlone=F,
#      width = 6, 
#      height = 8,
#      )
# 
# print(screegrid)
# 
# dev.off()

#Title and cap

```

# Balance tests

```{r Balance}

#X  = a + b*Treat +e

B1 <- lm(Samplenøgle %>% as.factor ~ age + gender + edu + income + beskæftigelse %>% as.factor + party + ideology + region, data = drpanel)

B1 <- lm(age ~ Samplenøgle, data = drpanel)
B2 <- glm(gender ~ Samplenøgle, data = drpanel, family = "binomial")
B1 <- 

stargazer(B2, type = "text")

drpanel$Samplenøgle <- drpanel$Samplenøgle %>% as.factor

drpanel$Samplenøgle <- relevel(drpanel$Samplenøgle, ref = "Ingen Treatment")
drpanel$Samplenøgle <-factor(drpanel$treatment, c("Ingen Treatment", "Kontrol",         "Geografisk","Risiko","Social","Tid", "Kontrol+håb", "Geografisk+håb","Risiko+håb", "Social+håb","tid+håb"))

test <- multinom(Samplenøgle ~ age + gender + eduyears + ideology, data = drpanel)
relratio <- exp(coef(test))

stargazer(test, 
          type="latex", 
          coef = list(relratio),
          star.char = c("\\dag", "*", "**", "***"),
          star.cutoffs = c(.1, .05, .01, .001),
          notes = c("† p<0,1; * p<0,05; ** p<0,01; *** p<0,001", "Referencegruppe er 'ingen stimulus'"),
          notes.append=F,
          p.auto = F,
          title = "Relative risikoratioer for at blive tildelt gruppe",
          header = F,
          single.row = F,
          intercept.bottom = F, 
          dep.var.caption  = "Stimuli",
          dep.var.labels.include = T,
          dep.var.labels = c("Kontrol", "Geografisk", "Risiko", "Social", "Tid","Kontrol", "Geografisk", "Risiko", "Social", "Tid"),
          covariate.labels = c("Konstant", "Alder, år", "Køn, mand", "Uddannelse, år", "Ideologi"),
          out = "DRmultinom.tex",
          label="DRMultinom",
          decimal.mark = ",",
          font.size = "small",
          column.sep.width = "-5pt",
          float.env = "sidewaystable",
          column.labels   = c("Intet håb", "Håb"),
          column.separate = c(5, 5))

z <- summary(test)$coefficients/summary(test)$standard.errors

p <- (1 - pnorm(abs(z), 0, 1)) * 2

```

#Running models

```{r test some assumptions}
resdf <- data.frame(matrix(ncol = 2, nrow = 1516))
resdfMTurk <- data.frame(matrix(ncol = 2, nrow = 524))

resdf$dr <- resid(m1_concern_index1)
resdfMTurk$MTurk<- resid(concern_m1)
qqnorm(res)
qqline(res)

plot(density(resid(concern_m1)))

##ggplot
ggqqDR <- ggplot(resdf, aes(sample = dr)) + stat_qq(alpha=0.05) + stat_qq_line() + theme_minimal() + labs(x="Teoretisk normalfordeling", y = "Residualer", title = "DR Panel")
ggqqMTurk <- ggplot(resdfMTurk, aes(sample = MTurk)) + stat_qq(alpha=0.05) + stat_qq_line() +theme_minimal() + labs(x="Teoretisk normalfordeling", y = "Residualer", title = "MTurk")

tikz(file = "qqplots.tex",
     standAlone=F,
     width = 6,
     height = 3,
     sanitize =T,
     )

print(plot_grid(ggqqDR, ggqqMTurk, 
          ncol =2))

dev.off()

ggdensDR <- ggplot(resdf, aes(x=dr)) + geom_density() + theme_minimal() + stat_function(fun = dnorm, n = 101, args = list(mean = 0, sd = .15)) + ylab("") + scale_y_continuous(breaks = NULL) + xlab("DR Panelet")

ggdensMTurk <- ggplot(resdfMTurk, aes(x=MTurk)) + geom_density() + theme_minimal() + stat_function(fun = dnorm, n = 101, args = list(mean = 0, sd = .15)) + ylab("") +   scale_y_continuous(breaks = NULL) + xlab("MTurk")

tikz(file = "densitetskurver.tex",
     standAlone=F,
     width = 6,
     height = 3,
     sanitize =T,
     )

print(plot_grid(ggdensDR, ggdensMTurk, 
          ncol =2))

dev.off()

```


```{r Do different videos impact concern?}
drpanel$male <- drpanel$gender %>% as.numeric()
drpanel$republican <- drpanel$right %>% as.numeric()
drpanel$Stimulus <- recode(drpanel$treatment5, Kontrol = 'Control', 
                                       Geografisk = 'Geographical',
                                       Risiko = 'Hypothetical',
                                       Tid = 'Time')

m1_concern_index1 <- lm(concern_index ~ Stimulus, data = drpanel)

HC_m1_concern_index1 <- coeftest(m1_concern_index1,vcov=vcovHC(m1_concern_index1))


m2_concern_index1 <- lm(concern_index ~ Stimulus + male + age + republican + eduyears, data = drpanel)
HC_m2_concern_index1 <- coeftest(m2_concern_index1,vcov=vcovHC(m2_concern_index1))

stargazer(m1_concern_index1,m2_concern_index1,concern_m1, concern_m2, type="latex",
          star.char = c("\\dag", "*", "**", "***"),
          star.cutoffs = c(.1, .05, .01, .001),
          notes = c("† p<0,1; * p<0,05; ** p<0,01; *** p<0,001", "Referencegruppe for afstande er kontrolgruppen", "Heterogenitetsrobuste standardfejl"),
          notes.append=FALSE,
          title = "Effekt af psykologiske afstande på bekymring for klimaforandringerne",
          header = F,
          single.row = T,
          se = list(HC_m1_concern_index1[,"Std. Error"],HC_m2_concern_index1[,"Std. Error"],HC_concern_m1[,"Std. Error"],HC_concern_m2[,"Std. Error"]),
          intercept.bottom = FALSE,
          dep.var.caption  = "Bekymring",
          dep.var.labels = "Bekymring",
          dep.var.labels.include = F,
          out = "bekymring.tex",
          label="bekymring",
          decimal.mark = ",",
          font.size = "footnotesize",
          column.sep.width = "-5pt", 
          covariate.labels = append(c("Konstant", "Geografisk", "Social", "Hypotetisk", "Tid", "Ingen stimuli", "Mand", "Alder, år", "Højreorienteret", "Uddannelse, år"),levels(MTurk$eduOrdered)[-1]),
          column.labels   = c("DR Panelet", "MTurk"),
          column.separate = c(2, 2))


m1_concern_index3 <- lm(concern_index3 ~ treatment5 + eduyears + treatment5*eduyears, data = drpanel)

m1_overall <- lm(Concern1 ~ treatment5, data = drpanel)
HC_m1_overall <- coeftest(m1_overall,vcov=vcovHC(m1_overall))


m1_social <- lm(Concern2_reverse ~ treatment5, data = drpanel)
#Concern2: Geografisk afstand øger fornemmelse af, at klimaet vil ramme vedkommende.
#Besynderligt! 7 procentpoint
HC_m1_social <- coeftest(m1_social,vcov=vcovHC(m1_social))


m1_geo <- lm(concern3 ~ treatment5, data = drpanel)
#Geography is again (marginally significant). 3,7 point
HC_m1_geo <- coeftest(m1_geo,vcov=vcovHC(m1_geo))


m1_uncertainty <- lm(Concern4_reverse ~ treatment5, data = drpanel)
#Time is marginally significant. Decreases by 5,2 points.
HC_m1_uncertainty <- coeftest(m1_uncertainty,vcov=vcovHC(m1_uncertainty))


m1_time <- lm(concern5 ~ treatment5, data = drpanel)
HC_m1_time <- coeftest(m1_time,vcov=vcovHC(m1_time))


stargazer(m1_geo,m1_social,m1_uncertainty,m1_time, type="latex", out = "DR_underdimensioner_bekymring.tex",
          star.char = c("\\dag", "*", "**", "***"),
          star.cutoffs = c(.1, .05, .01, .001),
          notes = c("† p<0,1; * p<0,05; ** p<0,01; *** p<0,001", "Referencegruppe for afstande er kontrolgruppen", "Heterogenitetsrobuste standardfejl"),
          notes.append=FALSE,
          title = "Effekt af psykologiske afstande på bekymringsmæssige underdimensioner",
          header = F,
          single.row = T,
          se = list(HC_m1_geo[,"Std. Error"],HC_m1_social[,"Std. Error"],HC_m1_uncertainty[,"Std. Error"],HC_m1_time[,"Std. Error"]),
          p=list(HC_m1_geo[,"Pr(>|t|)"],HC_m1_social[,"Pr(>|t|)"],HC_m1_uncertainty[,"Pr(>|t|)"],HC_m1_time[,"Pr(>|t|)"]),
          intercept.bottom = FALSE,
          dep.var.caption  = "Bekymringsmæssig underdimension",
          dep.var.labels = c("Geografisk", "Social", "Hypotetisk", "Tid"),
          dep.var.labels.include = T,
          covariate.labels = c("Konstant", "Geografisk", "Social", "Hypotetisk", "Tid", "Ingen treatment"),
          label="bekymringerDR",
          decimal.mark = ",",
          font.size = "footnotesize",
          column.sep.width = "-5pt")
```

```{r Do different videos impact efficacy?}

m1_efficacy_index <- lm(efficacy_index ~ Stimulus, data = drpanel)
HC_m1_efficacy_index <- coeftest(m1_efficacy_index,vcov=vcovHC(m1_efficacy_index))

m2_efficacy_index <- lm(efficacy_index ~ Stimulus + age + male + republican +eduyears  + income, data = drpanel)

HC_m2_efficacy_index <- coeftest(m2_efficacy_index,vcov=vcovHC(m2_efficacy_index))


m1_efficacy_index3 <- lm(drpanel$efficacy_index3 ~ treatment5 + age + gender + edu + right + income, data = drpanel)

stargazer(m1_efficacy_index, m2_efficacy_index, type="latex",
          star.char = c("\\dag", "*", "**", "***"),
          star.cutoffs = c(.1, .05, .01, .001),
          notes = c("† p<0,1; * p<0,05; ** p<0,01; *** p<0,001", "Referencegruppe for afstande er kontrolgruppen", "Heterogenitetsrobuste standardfejl"),
          se = list(HC_m1_efficacy_index[,"Std. Error"],HC_m2_efficacy_index[,"Std. Error"]),
          notes.append=FALSE,
          title = "Effekt af psykologiske afstande på følelse af efficacy",
          header = F,
          single.row = T,
          intercept.bottom = FALSE, 
          dep.var.caption  = "Efficacy",
          dep.var.labels = "Efficacy",
          dep.var.labels.include = F,
          covariate.labels = append(c("Konstant", "Geografisk", "Social", "Risiko", "Tid", "Ingen stimulus", "Alder", "Køn, mand", "Uddannelse, år", "Højreorienteret"),levels(drpanel$income)[-1]),
          out = "DRefficacyalder.tex",
          label="DRefficacyalder",
          decimal.mark = ",",
          font.size = "footnotesize",
          column.sep.width = "-5pt")

#There are no differences in either

```

```{r Do different videos impact intentions?}

m1_intention_index <- lm(intention_index ~ treatment5 + age + gender + edu + income + right, data = drpanel)

m1_efficacy_index3 <- lm(intention_index3 ~ treatment5 + age + gender + edu + income + right, data = drpanel)

#Nothing going on here

```

```{r Do the videos induce hope?}
mh1 <- lm(hope_index12 ~ treatment_hope, data = drpanel)
summary(mh1)

```

```{r Do the videos induce fear?}
mf1 <- lm(fear_index4 ~ treatment_hope, data = drpanel)
summary(mf1)

```

```{r Did the manipulations work?}
mc1 <- lm(check_faraway ~ treatment5, data = drpanel)
stargazer(mc1, type = "text")
#That worked

mc2 <- lm(check_likeme ~ treatment5, data = drpanel)
stargazer(mc2, type = "text")
#That worked too

mc3 <- lm(check_notlikely ~ treatment5, data = drpanel)
summary(mc3)
#That didn't work

mc4 <- lm(check_longtime ~ treatment5, data = drpanel)
summary(mc4)
#That didn't work either

mc5 <- lm(feels_hope ~ treatment_hope, data = drpanel)
stargazer(mc5, type = "text")
#That worked well

HC_mc1 <- coeftest(mc1,vcov=vcovHC(mc1))
HC_mc2 <- coeftest(mc2,vcov=vcovHC(mc2))
HC_mc3 <- coeftest(mc3,vcov=vcovHC(mc3))
HC_mc4 <- coeftest(mc4,vcov=vcovHC(mc4))


stargazer(mc1,mc2,mc3,mc4, type="latex", out = "manicheck.tex",
          star.char = c("\\dag", "*", "**", "***"),
          star.cutoffs = c(.1, .05, .01, .001),
          notes = c("† p<0,1; * p<0,05; ** p<0,01; *** p<0,001", "Referencegruppe for afstande er kontrolgruppen", "Heterogenitetsrobuste standardfejl"),
          notes.append=FALSE,
          title = "Manipulationstjek af intenderede effekter, DR Panelet",
          header = F,
          single.row = T,
          se = list(HC_mc1[,"Std. Error"],HC_mc2[,"Std. Error"],HC_mc3[,"Std. Error"],HC_mc4[,"Std. Error"]),
          intercept.bottom = FALSE,
          dep.var.caption  = "Overensstemelse mellem stimuli og oplevet afstand",
          dep.var.labels = c("Geografisk", "Social", "Hypotetisk", "Tid"),
          dep.var.labels.include = T,
          covariate.labels = c("Konstant", "Geografisk", "Social", "Hypotetisk", "Tid", "Ingen treatment"),
          label="manicheck",
          decimal.mark = ",",
          font.size = "footnotesize",
          column.sep.width = "-5pt")


```

You should make a figure showing this.

```{r Balance tables}

Balance1 <- drpanel %>%
  group_by(treatment) %>%
  summarize(AvgAge = mean(age, na.rm = T),
            SDAge = sd(age, na.rm = T),
            minAge = min(age, na.rm = T),
            maxAge = max(age, na.rm = T),
            
            PctWomen = mean(gender %>% as.numeric -1, na.rm = T),
            
            PctRight = mean(right, na.rm = T),
            sdRight = sd(right, na.rm = T),
            
            ) %>% cbind(drpanel %>% count(Samplenøgle))

#Count by income
inccount <- drpanel %>%
    group_by(income) %>% count(treatment)

inccount_wide <- spread(inccount, 
       key = treatment,
       value = n)

inccount_wide$income <-str_replace_na(inccount_wide$income, replacement = "NA")
inccount_wide <- inccount_wide %>% remove_rownames %>% column_to_rownames(var="income")

#Count by education
educount <- drpanel %>%
    group_by(edu) %>% count(treatment)

educount_wide <- spread(educount, 
       key = treatment,
       value = n)

educount_wide$edu <-str_replace_na(educount_wide$edu, replacement = "NA")
educount_wide <- educount_wide %>% remove_rownames %>% column_to_rownames(var="edu")

#Long to wide
spread(Balance1, 
       key = treatment,
       value = n)

balance_out_wide <- data.frame(Balance1 %>% t)

colnames(balance_out_wide) <- balance_out_wide[1,]
balance_out_wide <- balance_out_wide[-c(1),]

#Divider med totaler for procenter
educount_wide_pct <- educount_wide %>% dplyr::select(names(educount_wide)[1])*100/educount_wide %>% dplyr::select(names(educount_wide)[1]) %>% colSums()

for(i in 2:ncol(educount_wide)){
educount_wide_pct <- cbind(educount_wide_pct, educount_wide %>% dplyr::select(names(educount_wide)[i])*100/educount_wide %>% dplyr::select(names(educount_wide)[i]) %>% colSums())
    
}

educount_wide_pct <- rbind(educount_wide_pct, colSums(educount_wide_pct))

inccount_wide_pct <- inccount_wide %>% dplyr::select(names(inccount_wide)[1])*100/inccount_wide %>% dplyr::select(names(inccount_wide)[1]) %>% colSums()

for(i in 2:ncol(inccount_wide)){
inccount_wide_pct <- cbind(inccount_wide_pct, inccount_wide %>% dplyr::select(names(inccount_wide)[i])*100/inccount_wide %>% dplyr::select(names(inccount_wide)[i]) %>% colSums())

}

inccount_wide_pct <- rbind(inccount_wide_pct, colSums(inccount_wide_pct))

#End

tutti <- rbind(balance_out_wide,inccount_wide_pct,educount_wide_pct)
tutti[] <- lapply(tutti, as.numeric)

tutti %>% mutate(across(where(is.numeric), ~ round(., 2)))


#Lav balancetabel
print(xtable(tutti, caption="Balancetabel", digits=2, type = "latex"), file = "balance.tex")
kable(tutti,format = "latex", digits = 2, caption="Balancetabel")

```

```{r Excel as LaTeX}

strata <- read_excel("220009 - Treatmentfordeling og Samplestratifikation.xlsx", 
    sheet = "Stratifisering")

strata$`Fordling i befolkning %` <- strata$`Fordling i befolkning %`*100
strata$`Samplestørrelse %` <- strata$`Samplestørrelse %`*100
strata$Forskel <- strata$`Samplestørrelse %` - strata$`Fordling i befolkning %`


kable(strata %>% dplyr::select(c(1,3,5,6)),format = "latex", digits = 2)

```

```{r Rodeboks}

G1 <- drpanel %>%
  group_by(treatment5) %>%
  summarize(AvgConcern = mean(concern3, na.rm = T),
            SDConcern = sd(concern3, na.rm = T)) %>% 
  mutate(treatment5 = fct_reorder(treatment5, AvgConcern))


```

